Phase coupling estimation from multivariate phase statistics 
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Coupled oscillators are prevalent throughout the physical world. Dynamical system formulations of weakly 
coupled oscillator systems have proven effective at capturing the properties of real-world systems. However, 
these formulations usually deal with the 'forward problem': simulating a system from known coupling pa- 
rameters. Here we provide a solution to the 'inverse problem': determining the coupling parameters from 
measurements. Starting from the dynamic equations of a system of coupled phase oscillators, given by a non- 
linear Langevin equation, we derive the corresponding equilibrium distribution. This formulation leads us to 
the maximum entropy distribution that captures pair-wise phase relationships. To solve the inverse problem for 
this distribution, we derive a closed form solution for estimating the phase coupling parameters from observed 
phase statistics. Through simulations, we show that the algorithm performs well in high dimensions (d=100) 
and in cases with limited data (as few as 100 samples per dimension). Because the distribution serves as the 
unique maximum entropy solution for pairwise phase statistics, the distribution and estimation technique can be 
broadly applied to phase coupling estimation in any system of phase oscillators. 

PACS numbers: 05.45.Xt, 05.10.Gg, 87.19.ln, 89.75.-k 



Many complex natural phenomena can be modeled as net- 
works of coupled oscillators. Examples can be drawn from 
the physical, chemical, and biological world. Oscillator mod- 
els have been effective at describing the dynamics of coupled 
pendula, coupled Josephson junctions, reaction diffusion sys- 
tems, circadian rhythms, oscillating neural networks, and even 
the coupling of firefly luminescence (see e.g. IIIlEllSllllO). 

In many systems, coupling topology and the strength of in- 
teraction between network elements is of central scientific in- 
terest. However, network coupling often can not be measured 
directly and must be inferred from measurements. Therefore 
the inverse problem, or inferring network coupling from mea- 
surements, is of central importance. 

In statistical mechanics, the inverse problem is typically 
solved by proposing a probability distribution and estimating 
the distribution's parameters from measurements. A natural 
choice for the estimation, a highly under-determined prob- 
lem, is given by the unique maximum entropy distribution that 
reproduces the statistics of the measurements |6|. A num- 
ber of such distributions and estimation techniques are used 
throughout the science and engineering communities. In the 
real- valued case the multivariate Gaussian distribution, and in 
the binary case the Ising model, serve as widely used mul- 
tivariate maximum entropy distributions consistent with sec- 
ond order statistics. Each of these cases has well known es- 
timation techniques for inferring the distribution's parameters 
from observations. The availability of these techniques has led 
to a number of applications, e.g. the Ising model and its cor- 
responding estimation techniques have been used to infer the 
coupling in networks of retinal ganglion cells liZlEl- However, 
for the phase variables that are of interest in networks of oscil- 
lators there has been little work on providing a corresponding 
multivariate probabilistic distribution, or deriving estimation 
techniques to infer the distribution's parameters from data. 



In this Letter, we provide a solution to the inverse prob- 
lem for systems of coupled phase oscillators. We begin by 
presenting the Langevin dynamics for a generalized form of 
the Kuramoto model of coupled phase oscillators. Solving for 
the equilibrium distribution yields a multivariate probability 
distribution of coupled phase variables. This probabilistic for- 
malism allows us to derive a novel estimation technique for 
the coupling terms from phase variable measurements. We 
show that this technique performs robustly with limited data 
and in high dimensions. 

Consider a network of d identical coupled oscillators with 
intrinsic frequency co. In the limit of weak coupling, the am- 
plitude of the oscillators can be assumed to be constant and 
the equations of motion can be formulated in terms of d phase 
variables < 6i < 27r, i = 1, . . . , (i. A popular choice for 
the dynamics of such a system is given by the Kuramoto 
model 1 2 1, which has constant coupling between oscillators. 
We can generalize this model to include inhomogeneous cou- 
plings and inhomogeneous phase offsets between oscillators. 
The dynamic equation is then given by 

^ 

-Q-^0,{t)=uo - ^sm{Oi{t)-Oj{t)-^ij) + Ti{t), (1) 

where f^ij IS the coupling strength and /lij is the preferred 
phase between two oscillators i and j. We only consider the 
case of symmetric coupling (hZij = Kji, fiij = —fiji). The noise 
fluctuations, Ti{t), are zero mean Gaussian distributed with S 
covariance functions and variance corresponding to the 
temperature of the system: 

{r,{t)) = 0, (r,(t)r,(tO) = 2p-'5,,5{t - 1') . (2) 

The equations of motion ([T]) for our system of coupled os- 
cillators can be considered as a nonlinear Langevin equation 
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describing Brownian motion on a (i-torus in the presence of 
the potential E{6) given by 
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where 6 is now a d-dimensional vector with components Oi. 
Note that by applying the transformation Oi {t) = Oi (t) — cot 
to ([T]) we can assume cc; = without loss of generality. 

By changing the coordinates from the angular representa- 
tion, 6, to the complex representation, {x G | = 1} 
with components Xk = e*^^ , we can rewrite eq. ^ more com- 
pactly as the (real- valued) quadratic Hermitian form: 

£^(x) = -^x^Kx, (4) 

where K is a d x d Hermitian matrix with elements Kji^ = 
Kjkc'^^^'^. This energy function ^ is closely related to the 
XY-model, which only has homogeneous nearest neighbor 
couplings (kij = const.) and no phase offsets (/i^^=0). This 
generalization is analogous to the extension of the homoge- 
neous Ising model to spin glasses. 

It is known (see e.g. [9l) that the probability density 
p{6{t)^t) of a system governed by Langevin dynamics 
evolves according to the Fokker-Planck equation 
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with drift and diffusion coefficients given by 
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Since the drift coefficient Di is a gradient field and the dif- 
fusion coefficient Dij is constant, we can solve the Fokicer- 
Planck equation ^ for the stationary solution in closed form 
and obtain 



Z{K) 



(7) 



with the the energy function E{0) given by ([s]), and partition 
function Z(K). 

We wish to solve the inverse problem for the general case 
of coupled phase oscillators in equation (|7]). Stated explicitly, 
the problem is to infer the distribution's parameters (coupling 
terms K,ij and phase offsets jj^ij) from measurements of the 
network's state, 6. 

The inverse problem is typically solved by following a max- 
imum likelihood estimation procedure. Given the likelihood 
function q{6) and the data distribution p{6), the maximum 
likelihood of the observed data with respect to the distribution 
parameters can be computed by setting the derivative of the 
log-likelihood function to zero. 
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= 0, (8) 



where (. . denotes the expectation value taken over the 
distribution q{6). In our situation, a closed form solution to 
this equation does not exist. However, we can find a solu- 
tion by iteratively descending the gradient. This procedure 
has a number of drawbacks: the procedure is inherently iter- 
ative, estimating the expectation under the estimated distribu- 
tion q{0) in equation ^ involves a computationally expensive 
sampling procedure, and the sampling procedure may suffer 
from a variety of problems due to the landscape of the energy 
function. 

To avoid the pitfalls of the maximum likelihood approach, 
we now derive a closed form solution to the inverse prob- 
lem for phase coupled oscillators. We use the score-matching 
method introduced by Hyvarinen ITOl HB. Score-matching 
allows the fitting of probability distributions of the exponen- 
tial form for real-valued data without computation of the nor- 
malization constant Z. If the energy depends linearly on 
the distribution parameters, the solution can be computed in 
closed form by setting the derivative of the score function to 
zero fTT\. 

We follow this approach to estimate the distribution param- 
eters for our distribution in equation ([7]). The score matching 
estimator of K is given by K = argminx Jsm(K) and the 
score function Jsm(K) is given by 



Jsm(K) = (^[V^^(^)][V^^(^)]^ 



with the expectation value, (...), taken over the data distribu- 
tion. Using the quadratic form of the energy in ^ and the 
Jacobian Dij := dxi/dOj, we compute 
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The estimator K is computed by setting the derivative of the 
score function d/dKij Jsm(K) to zero. This produces a sys- 
tem of linear equations. 
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{^jiCik + SikCij - SjkCiijk — SilCiijk) Kki = 4 Ci^ 



(9) 

where the expectation values are defined as Cij = {xixD and 
Cijki = {xiXjX^Xi). We can solve this system of linear equa- 
tions using standard techniques. 

In the following, we show that phase coupling estimation 
recovers the parameters of simulated Kuramoto models. We 
begin by simulating a system of four oscillators using equa- 
tion ([T]) with couplings shown in Fig. [T^. Given samples of 
the simulated phase variables 6, we compute the correlations, 
Cij, and required four-point functions, Cijki, and invert the 
linear system in equation ([9]). This produces an estimate of 
the coupling parameters. Phase coupling estimation recovers 
the true coupling parameters as shown in Fig.[T]D. 
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FIG. 1 : Phase coupling estimation recovers the true coupHng from 
measurements, (a) Diagram of four coupled phase oscillators. 

(b) Estimated coupling parameters using phase coupling estimation. 

(c) Pair- wise phase correlations of the oscillators in a. Phase corre- 
lations are parameterized by the mean direction vector rue^^^^ — 
(e*^^e~*^^) with amplitude Vki and angle /\ki- (d) Empirical phase 
distribution of the phase difference between oscillators 1 and 2. The 
empirical distribution is highly concentrated in the difference of the 
phases, while the marginals are flat (not shown). Phase correlations 
and estimated couplings are given, (e) Distribution of phase differ- 
ences of oscillators 1 and 4, which are not directly coupled, (f) Uni- 
form distribution of phase differences of oscillators 2 and 3, which 
are directly coupled. 



We would like to point out that the pairwise statistics in 
systems of coupled phase oscillators are only indirectly re- 
lated to the coupling parameters. Phase correlations, a pair- 
wise measurement often used to characterize oscillator sys- 
tems, have a direct relationship to the marginal distribution of 
phase differences but not to coupling parameters. The form of 
the marginal distribution can be derived by examination of the 
individual factors in the definition 



p{ek-ei) 



n 

i.3 = ^ 



exp [nij cos{Oi — Oj — i^ij)] dO' 



d-2 



(10) 

with i k^l. 
+ Oi , all terms 

in ([To]) either depend on the phase difference O^—Oi, or are 



in which the integration is over all phases 0^ 
After applying the variable substitution Oi = Oi 



independent of 0^ and Oi. The independent terms integrate to 
a constant and the remaining terms combine to a von Mises 
distribution in the pairwise phase difference 



p{0k-0i) 
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^jki cos{0k-0i-Aki) 
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with mean phase Aki and concentration parameter jki- The 
concentration parameter jki can be obtained by numerically 
solving the equation rki = h (jki) / ilki) and the normaliza- 
tion constant Z{'yki) is given by Z{jki) = 2 7r/o(7/cO- ^o(^) 
and Ii{x) denote the modified Bessel functions of zeroth and 



first order, respectively. The value of jki is related to the cou- 
pling parameters K through equation ([TO]). Therefore, there is 
a non-trivial relationship between the measured phase corre- 
lations and the coupling parameters. 

Because of this non-trivial relationship, pair- wise measure- 
ments can often lead to false interpretations of the true cou- 
pling. We show the measured phase correlations of our four 
oscillator system in Fig.[T]:. Phases Oi and O4 show clear cor- 
relations even though they are not directly coupled to each 
other (see Fig.[T^). Conversely, phases O2 and Os are strongly 
coupled but are uncorrected (see Fig.[T]r). 

We now systematically analyze the performance of phase 
coupling estimation: the ability of the technique to recover 
the distribution parameters from data. The procedure is as fol- 
lows. We begin by sampling a set of distribution parameters 
K. Given these parameters we then sample phase variables 
by numerically integrating ([T]). We then estimate the parame- 
ters given the sampled data using equation ([9]). The real and 
imaginary entries of the complex matrix K are sampled from 
a normal distribution: Re{Kij}^ lm{Kij} ~ A/'(0, 1) and the 
diagonal entries are set to zero: Ka = 0. Note that this pro- 
duces a dense coupling matrix. 
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FIG. 2: Phase coupling estimation for a system of 16 coupled os- 
cillators with random coupling, (a) True coupling matrix K: true 
system parameters for d = 16 (first row, element-wise amplitude; 
second row, element- wise angle with alpha channel scaled by the am- 
plitude of the corresponding element, best viewed in color), (b) es- 
timated coupling matrix K: estimated parameters recovered from 
2560 time samples of 6 using equation (c) estimation error: first 
row element-wise mse (note scaling), second row estimation error 
measurements, mse and Q.95 (see text for definition). 

In the first column of Fig. [2j we graphically display the 
element-wise amplitude and phase of a sample matrix K 
where d=16. Using this matrix we sampled N = 2560 phase 
vectors. The recovered parameters are shown in the second 
column of Fig. [2] While it is clear that these matrices are 
visually similar, we quantified the error using two different 
metrics. First we calculate the mean-squared-error of the ma- 
trix elements, mse = j where Kij is the 
estimated parameter. In the third column of Fig. [2] we display 
the element- wise error before averaging. We also computed a 
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metric indicating the quality of the recovered parameters bor- 
rowed from Ref. 1 12|: (5.95 = ^ j u{l— .95AKij), where 

AKij = \Kij — Kij\/2K^s,x, u is the Heaviside step func- 
tion, and i^Trnax is the maximum absolute value of all matrix 
entries Kij and Kij. For the example in Fig. [2} mse = 0.02, 
and Q.95 = 0.75. 

We computed these error metrics over a range of dimen- 
sions and samples per dimension. The error metrics for each 
dimension and samples per dimension were averaged over 20 
trials and are plotted in Fig. [3] The algorithm achieves highly 
accurate parameter recovery for as few as 100 samples per di- 
mension and achieves full recovery of parameters as the num- 
ber of samples per dimension reaches 1000. This indicates 
that recovery of true parameters is quite feasible in many real 
world settings. 
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FIG. 3: Performance of phase coupling estimation, (a) mean- 
squared-error, mse, metric as a function of samples per dimension 
for various dimensions, d = 8, 16, 32, 64, 100. (b) Q.95 metric. The 
example displayed in Fig.|2]is indicated by the black diamond. Val- 
ues are averaged over 20 trials. For visual clarity, standard error bars 
are only displayed for d = 8 and d = 100. Phase coupling estima- 
tion accurately recovers the system parameters with only 100 sam- 
ples per dimension and achieves nearly perfect recovery with 1000 
samples per dimension. 

In this Letter we have introduced a closed form solution to 
the inverse problem for systems of coupled phase oscillators 
using a maximum entropy approach. We close by pointing out 
the relation of our work to other attempts at solving the inverse 
problem for coupled phase oscillators. Previous examples of 
probabilistic distributions of phase variables only characterize 
univariate or bivariate distributions (see e.g. (13] [141). Distri- 
butions similar to equation ^ have not been extended to di- 
mensions beyond d = 2 [15^ . Common multivariate phase 



distributions do not capture the statistics that are relevant for 
coupled phase oscillators. Most notably the von Mises-Fisher 
distribution only captures a unimodal phase distribution on the 
hyper-sphere, which does not produce unimodal concentra- 
tions in the differences of phase variables. One of the most 
relevant proposals is the estimation procedure of Timme lfT2l . 
While the method of Ref. 1 12] successfully recovers the cou- 
pling parameters, it requires repeated intervention, which may 
not be feasible in many real- world experiments. 

Phase coupling estimation can potentially provide a contri- 
bution in a variety of situations of scientific interest. Because 
our phase coupled estimation technique derives the unique 
maximum entropy solution, it serves as the least biased esti- 
mate of the system possible and can be used when the true dy- 
namics of the system are unknown. Such situations are preva- 
lent in neuroscience where oscillations are thought to mitigate 
cognitive processes but the form of the underlying dynamical 
system is unknown. This field has lacked a suitable proce- 
dure for estimating phase interactions and has largely relied 
on pair-wise statistical measurements (see e.g. 1 17 1). 
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